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Abstract 

We present a simple and powerful method for extrapolating finite- volume 
Monte Carlo data to infinite volume, based on finite-size-scaling theory. We 
discuss carefully its systematic and statistical errors, and we illustrate it using 
three examples: the two-dimensional three-state Potts antiferromagnet on the 
square lattice, and the two-dimensional 0(3) and O(oo) u-models. In favorable 
cases it is possible to obtain reliable extrapolations (errors of a few percent) 
even when the correlation length is 1000 times larger than the lattice. 
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No matter how powerful computers get, physicists will always want to study 
problems that are too difficult for the computers at hand. For example, in statis- 
tical mechanics and quantum field theory, physicists want to push to ever larger 
correlation lengths £. But Monte Carlo simulations must perforce be carried out on 
lattices of finite linear size L (limited by computer memory and speed); the data are 
then extrapolated to the infinite- volume limit L = oo. Obviously this extrapolation 
— which is based on the theory of finite-size scaling (FSS) [1,2,3] — is feasible in 
practice only if £/L is not too large. But how large? 

In this Letter we present a simple and powerful method for performing the 
extrapolation to L = oo, and discuss carefully its systematic and statistical errors. 
We illustrate the method using three examples: the two-dimensional three-state 
Potts antiferromagnet on the square lattice [4], and the two-dimensional 0(3) and 
0(oo) cr-models [5,6]. We have found — much to our surprise — that in favorable 
cases it is possible to obtain reliable extrapolations (errors of a few percent) at £/L 
as large as 10-1000. 

Consider, for starters, a model controlled by a renormalization-group (RG) fixed 
point having one relevant operator. Let us work on a periodic lattice of linear size 
L. Let £(f3,L) be a suitably defined finite- volume correlation length [7], and let O 
be any long-distance observable (e.g. the correlation length or the susceptibility). 
Then finite-size-scaling theory [1,2,3] predicts that [9] 

° { ^ L) foU(P,oo)/L) + 0(C u ,L-«) , (1) 



where fo is a universal function and u is a correction-to-scaling exponent. It follows 
that if s is any fixed scale factor (usually we take s = 2), then 

= F (t(P,L)/L) + 0(C,L-») , (2) 

where Fq can be expressed in terms of fo and f^. 

Our method proceeds as follows [10]: Make Monte Carlo runs at numerous pairs 
(/3,L) and (/3 } sL). Plot 0(/3 } sL)/ 0(/3 } L) versus £(f3, L)/ L, using those points 
satisfying both £(f3,L) > some value ^ m in and L > some value L m i n . If all these 
points fall with good accuracy on a single curve — thus verifying the Ansatz (2) 
for £ > £ m in, L > L m i n — choose a smooth fitting function Fq- Then, using the 
functions F% and Fq, extrapolate the pair (£, O) successively from L — > sL — > s 2 L — > 
oo. 

We have chosen to use functions Fq of the form 

Fo(x) = 1 + a x e~ x l x + a 2 e~ 2 l x + ... + a n e~ n l x . (3) 

This form is partially motivated by theory, which tells us that F(x) — > 1 exponen- 
tially fast as x — > [14]. Typically a fit of order 3 < n < 9 is sufficient; we increase n 
until the % 2 of the fit becomes essentially constant. The resulting % 2 value provides 
a check on the systematic errors arising from corrections to scaling and/or from the 
inadequacies of the form (3). 
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The statistical error on the extrapolated value of 0qo(/3) = 0(/3 } oo) comes from 
three sources: 

(i) Error on 0(/3 } L) } which gets multiplicatively propagated to Cqo- 

(ii) Error on which affects the argument x = £(f3,L)/L of the scaling 
functions F% and Fq- 

(iii) Statistical error in our estimate of the coefficients ai, . . . , a n in F% and Fq- 

The errors of type (i) and (ii) depend on the statistics available at the single point 
(/3, L), while the error of type (iii) depends on the statistics in the whole set of runs. 
Errors (i)+(ii) [resp. (i)+(ii)+(iii)] can be quantified by performing a Monte Carlo 
experiment in which the input data at (/3, L) [resp. the whole set of input data] are 
varied randomly within their error bars and then extrapolated [15,16,17]. 

The discrepancies between the extrapolated values from different lattice sizes at 
the same f3 — to the extent that these exceed the estimated statistical errors — 
indicate the presence of systematic errors and thus the necessity of increasing L m i n 
and/or ^ m in and/or n. 

A figure of (de)merit of the method is the relative variance on the extrapolated 
value (9oo(/3), multiplied by the computer time needed to obtain it [18]. We expect 
this relative variance-time product [for errors (i)+(ii) only] to scale as 

RVTP(/3,L) « iU0 +z > nt ' o Go(ioo(P)/L) , (4) 

where d is the spatial dimension and Zi nt ,o 1S the dynamic critical exponent of the 
Monte Carlo algorithm being used; here Go is a combination of several static and 
dynamic finite-size-scaling functions, and depends both on the observable O and 
on the algorithm but not on the scale factor s. As £,co/L tends to zero, we expect 
Go to diverge as (^ 00 /F)~ d (it is wasteful to use a lattice L ^> £oo)- As ^oo/L tends 
to infinity, we expect Go ~ (Cco/L) p [19], but the power p can be either positive or 
negative. If p > 0, there is an optimum value of ^/L; this determines the best 
lattice size at which to perform runs for a given f3. If p < 0, it is most efficient to use 
the smallest lattice size for which the corrections to scaling are negligible compared 
to the statistical errors. [This neglects errors of type (iii); the optimization becomes 
much more complicated if they are included.] 

Our first example [4] is the two-dimensional three-state Potts antiferromagnet 
on the square lattice, which is believed to have a critical point at f3 = oo [20]. We 
used the Wang-Swendsen-Kotecky cluster algorithm [21], which appears to have 
no critical slowing-down {r int ^ M 2 < 5 uniformly in f3 and L) [4]. We ran on 
lattices L = 32, 64, 128, 256, 512, 1024, 1536 at 153 different pairs (/3, L) in the range 
5 ^ iJjoo ^ 20000. Each run was between 2 X 10 5 and 2.2 X 10 7 iterations, and the total 
CPU time was modest by our standards (about 2 years on an IBM RS-6000/370). 
We took ^min = 10 and L m i n = 64 and used a quintic fit in (3); the result for F% is 
shown in Figure 1 (% 2 = 56.47, 84 DF, level = 99%). The extrapolated values from 
different lattice sizes at the same f3 agree within the estimated statistical errors 
(% 2 = 35.79, 93 DF, level > 99%): see Table 1 for an example. The result for G(_ is 
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shown in Figure 2; the errors are roughly constant for £,co/L > 0.4 but rise sharply 
for smaller ^/L. The theoretical exponent p computed [19] from the fitted F% is 
equal to p = —0.09 ± 0.06; the curve suggests that p = 0. In practice we were able 
to obtain ^ to an accuracy of about 1% (resp. 2%, 3%, 5%) at ^ 1000 (resp. 
2000, 5000, 10000). 

Next let us consider [5,6] the two-dimensional 0(3) cr-model. We used the 
Wolff embedding algorithm with standard Swendsen-Wang updates [22,23,8]; again 
critical slowing-down appears to be completely eliminated. We ran on lattices 
L = 32,48,64,96,128,192,256,384,512 at 169 different pairs (J3,L) in the range 
20 ^ ijjoo ^ 10 5 . Each run was between 10 5 and 5 X 10 6 iterations, and the total 
CPU time was 7 years on an IBM RS-6000/370. We took ^ m in = 20 and used an 
eighth-order fit. There appear to be weak corrections to scaling (of order < 1.5%) 
in the region 0.3 ^ Cl/L ^ 0.7 for lattices with L < 64-96. We therefore chose 
conservatively L m i n = 128 for ^l/L < 0.7, and L m i n = 96 for ^l/L > 0.7. The 
result for F% is shown in Figure 3 (x 2 = 38.62, 62 DF, level = 99%). The result for 
G(_ is shown in Figure 4; at large £,co/L it decreases sharply, with a power p ~ — 2 
in agreement with theory [19]. In practice we obtained ^ to an accuracy of about 
0.8% (resp. 1.4%, 2.1%) at ^ ~ 10 3 (resp. 10 4 , 10 5 ). 

We also carried out a "simulated Monte Carlo" experiment for the O(N) cr- 
model at N = oo, by generating data from the exact finite- volume solution plus 
random noise of 0.1% for L = 64,96,128, 0.2% for L = 192,256 and 0.5% for 
L = 384,512 [which is the order of magnitude we attain in practice for 0(3)]. We 
considered 35 values of f3 in the range 20 < Coo £ 10 6 . We used £ min = 20 and 
Lmin = 64 (in fact much smaller values could have been used, as corrections to 
scaling are here very small) and a ninth-order fit; for two different data sets we get 
X 2 = 114 (resp. 118) with 166 DF. In practice we obtain ^ with an accuracy of 
0.6% (resp. 1.2%, 2%, 3%) at ^ w 10 3 (resp. 10 4 , 10 5 , 10 6 ). Here we can also 
compare the extrapolated values £^, tr (f3) with the exact values i e ^ act {fi) ■ Defining 

7L = E[Co*(/3) - CT ct (/3)]7^ 2 (/3), we find for the two data sets Tl = 17.19 (resp. 

P 

25.81) with 35 DF. Only 6 (resp. 9) points differ from the exact value more than 
one standard deviation, and none by more than two. 

Details on all of these models will be reported separately [4,6]. 

The method is easily generalized to a model controlled by an RG fixed point 
having k relevant operators. It suffices to choose k — 1 dimensionless ratios of 
long-distance observables, call them R = (i?i, . . . , Rk-i); then the function Fq will 
depend parametrically on R(f3 } L). In practice one can divide i?-space into "slices" 
within which Fq is empirically constant within error bars, and perform the fit (3) 
within each slice. We have used this approach to study the mixed isovector/isotensor 
cr-model, taking R to be the ratio of isovector to isotensor correlation length [5,6]. 

The method can also be applied to extrapolate the exponential correlation length 
(inverse mass gap). For this purpose one must work in a system of size L d ~ x X T 
with T ^> ^ exp (f3 } L) (compare [11]). 

We wish to thank Martin Hasenbusch and especially Jae-Kwon Kim for sharing 
their data with us, and for challenging us to push to ever larger values of £/L. The 
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e(L)/L 

Figure 1: 2L)/£(/?, L) versus £(f3,L)/L for the two-dimensional three-state 

Potts antiferromagnet. Symbols indicate L = 32 (+), 64 (x), 128 (□), 256 (O), 
512 (o)- Error bars are one standard deviation. Curve is a quintic fit in (3), with 
train = 10 and L m i n = 64. 
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Figure 2: Relative variance-time product [for errors (i)+(ii) only] divided by ^oo(/3) 2 , 
plotted versus ^ 00 (/3)/L, for two-dimensional three-state Potts antiferromagnet. 
Symbols indicate L = 32 (+), 64 (x), 128 (□), 256 (O), 512 (o), 1024 (*), 1536(*). 
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Figure 3: 2L)/£(/?, L) versus £(f3,L)/L for the two-dimensional 0(3) cr-model. 
Symbols indicate L = 32 (+), 48 (*), 64 (x), 96 (&), 128 (□), 192 ()^(), 256 (O). 
Error bars are one standard deviation. Curve is an eighth-order fit in (3), with 
Z min = 20 and L mm = 128 (resp. 96) for i{L)/L < 0.7 (resp. > 0.7). 
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0.1 1 10 100 1000 

Figure 4: Relative variance-time product [for errors (i)+(ii) only] divided by ^oo(/3) 2 , 
plotted versus ^ 00 (/3)/L, for two-dimensional 0(3) cr-model. Symbols indicate L = 
32 (+), 48 (*), 64 (x), 96 (&), 128 (□), 192 (ft), 256 (O), 384 (<>-), 512 (o). For 
comparison, the line shows the theoretical limiting slope —2. 
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Raw Data 


Extrapolated 


L 


Iterations 






oo) 


32 


2.2 x 10 7 


19.02 (0.01) 


90.51 ( 


1.01) 


64 


10 6 


35.52 (0.04) 


92.66 ( 


0.79) 


128 


10 6 


60.61 (0.09) 


93.17 ( 


0.42) 


256 


10 6 


84.69 (0.15) 


93.19 ( 


0.30) 


512 


5 x 10 5 


92.48 (0.33) 


92.89 ( 


0.38) 


1024 


2 x 10 5 


93.78 (1.17) 


93.78 ( 


1.16) 


mean 






93.13 ( 


0.26) 


(L > 64) 




X 2 = 0.85 (4 DF, level = 


= 93%) 



Table 1: Raw and extrapolated correlation lengths for the two-dimensional three- 
state Potts antiferromagnet at f3 = 3.5. Extrapolation based on ^ m in = 10 and 
Lmin = 64 and a quintic fit. For each extrapolated value we have reported the 
standard deviation of the estimate, including errors of all three types. The mean 
value and the % 2 have been computed taking into account the full covariance matrix 
[17]. 



11 



